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In a recent letter [Europhys. Lett. 97, 34002 (2012)], random matrix theory is introduced for long- 
range acoustic propagation in the ocean. The theory is expressed in terms of unitary propagation 
matrices that represent the scattering between acoustic modes due to sound speed fluctuations 
induced by the ocean's internal waves. The scattering exhibits a power-law decay as a function of 
the differences in mode numbers thereby generating a power-law, banded, random unitary matrix 
ensemble. This work gives a more complete account of that approach and extends the methods to 
the construction of an ensemble of acoustic timefronts. The result is a very efficient method for 
studying the statistical properties of timefronts at various propagation ranges that agrees well with 
propagation based on the parabolic equation. It helps identify which information about the ocean 
environment can be deduced from the timefronts and how to connect features of the data to that 
environmental information. It also makes direct connections to methods used in other disordered 
wave guide contexts where the use of random matrix theory has a multi-decade history. 



PACS numbers: 43.30.Ft, 43.30.Bp, 05.45.Mt 

I. INTRODUCTION 

Our interest in this paper is in long-range acoustic 
propagation that is vertically confined by the deep-ocean 
sound channel. In addition to the vertical wave guide 
confinement, there is important multiple scattering of 
the acoustic waves that is induced by local fluctuations 
in the ocean environment?^. A number of previous ex- 
periments were designed with an acoustic source placed 
centrally in the wave guide, i.e. near the so-called sound 
channel axis, and receivers placed in vertical arrays^'. 
The acoustic waves arrive at the receivers in timefronts 
that feature early arriving isolated branches that give 
way to complicated finales with the majority of the en- 
ergy. The timefronts have both structural and statistical 
properties which are the focus of this investigation. 

An improved understanding of the acoustic timefronts 
would lead to the possibility of deducing new informa- 
tion about the state of the ocean. However, the multiple 
scattering greatly complicates the analysis of experimen- 
tal data, and in particular is responsible for a partial loss 
of information about the ocean"™. It is still not fully 
resolved as to what information can be deduced from 
long-range propagation and how to extract it. 

Since the 1970's, numerical studies have often utilized 
paraxial optical approximations^ coupled with Fourier 
transform based projections to simulate acoustic time- 
fronts. It has advanced over the years into an indis- 
pensable computational tool, but it is not well suited 
to extracting analytic relationships between properties of 
the ocean and timefronts or providing insights that come 
readily with ray methods. It furthermore is still compu- 
tationally intensive enough that performing larger scale 
ensembles of calculations for statistical studies would be 
prohibitive. 

Although, the ocean wave guide has its own unique fea- 
tures, such as the timefront itself, it also possesses com- 
mon characteristics with many other wave guide prob- 
lems. Indeed, long-range ocean acoustics has long been 



recognized as a context in which quantum/ wave chaos 
is presenlPHUD Therefore it has strong parallels with 
physical systems such as mesoscopic disordered and bal- 
listic electronic conductor a 14 ! 15 ! and quantum point con- 
tacts^. For example, it was noted that quantum point 
contacts channel wave energy, called coherent branching 
in that context (not to be confused with branches of a 
timefront), through the behavior of its underlying rays^l 
in much the same way as occurs in long-range ocean 
acoustics^. Curiously, one of the key theoretical tools 
of quantum chaos, random matrix theory, was only re- 
cently initiated by us in the ocean acoustic propagation 
contextP^o] Partly this is due to the fact that tradi- 
tional random matrix theory does not straightforwardly 
translate there. In other acoustic contexts, e.g. certain 
problems in elasto-dynamics 2 ^, there exists a more di- 
rect application i n whi ch the featureless ensembles of 
Wigner and DysorP^^l apply without modification. In 
ocean acoustics, random matrix theory enters in terms 
of how amplitudes in the acoustic modes mix with prop- 
agation in range, but some dynamical information must 
be incorporated that is absent in the standard ensem- 
bles. The mixing is captured by the unitary propagation 
matrix of complex probability amplitudes for modal tran- 

sitiorPna. 

Random matrix theory is connected to a maximum 
entropy hypothesis 28 and is well suited for investigating 
what information can be deduced and how it appears in 
long-range acoustic propagation. Thus, there is a multi- 
fold purpose to this paper; i) to give a more complete 
description of how random matrix theory can be intro- 
duced into long-range acoustic propagation; ii) apply it 
for the first time to the full construction of timefronts; 
iii) investigate possible new connections between environ- 
mental information and experimentally accessible quan- 
tities; and iv) examine some of the statistical properties 
that emerge from this approach. 

This paper is organized as follows. The next section 
gives basic background material both for the ocean acous- 
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tics system and random matrix theory. This serves to 
provide the essential physical context, introduce possibly 
unfamiliar concepts from random matrix theory, spec- 
ify what will be taken into account, and define notation. 



Section III gives the derivation of the random matrix en- 
sembles to be used to model the acoustical system and 
points out which environmental information is accounted 
for in this first application and construction. Section [IV 
shows the comparison between timefronts produced with 
parabolic equation calculations and random matrix en- 
sembles. Both average and fluctuating properties are 
considered. We finish with a discussion of the merits 
of the random matrix approach and future possibilities. 



II. BACKGROUND - OCEAN ACOUSTICS AND 
RANDOM MATRIX THEORY 

The following four subsections specify the features of 
the long-range acoustics system which are taken into ac- 
count in this study. For our purposes, the idea is not to 
have the greatest realism and most accurate modeling of 
the ocean, but rather to have enough complexity to illus- 
trate how random matrix theory can be applied and be 
helpful. The final subsection contains some elementary 
background in random matrix theory since this may be 
unfamiliar to many readers. We summarize a few use- 
ful concepts such as ergodicity, statistical measures, and 
information. 



A. Parabolic equation 

The parabolic equation for the propagation of an 
acoustic wavefield ^(z, r; k) as a continuous wave^ with 
wavenumber /c, depth z, and range r through a sound 
speed potential V is 

~*(z, r;k) = - J^*(*> r; k) + V(z, r)tt(z, r; k) 

(1) 

This equation assumes that the propagation is essentially 
forward and thus the energy moves within a narrow band 
of angles relative to the sound channel axis and that there 
is no significant backscattering. 

The sound speed can be decomposed into the reference 
sound speed, cq, and fluctuations, 5c, about the reference: 
c(z, r) = co + Sc(z, r) with 5c(z, r) « Co, so the potential 
is approximated as follows: 



Sc(z, r) 




V(z,r) = -ll 



Our calculations use cq = 1.49 km/sec. 

B. Sound speed model 

The model potential V(z, r) takes the form 

V {z,r) = = V (z) + eV 1 (z,r) 

Co 



(2) 



(3) 



where Vo represents the relative change of the sound 
speed due to the vertical waveguide confinement and eV\ 
represents the relative sound speed fluctuations due to in- 
ternal waves, which are small perturbations on the wave 
guide 7 ; e is of order 10 -3 . In fact, other perturbations 
could easily be accounted for and added into eVi(z,r), 
but for the purposes of this work accounting for the in- 
ternal waves is sufficient and the overall physical context 
no less general. 

We take the vertical confinement to be expressed by 
Munk's canonical model 29 



V (z) 



B-y 



■rj(z) 



(4) 



where rj(z) = 2[z — z a ]/B, z a = 1 km is the sound channel 
axis, B is the thermocline depth scale giving the approx- 
imate width of the sound channel, and 7 is a constant 
representing the overall strength of the confinement. 

The fluctuations due to internal waves are captured 
with the statistical ensemble model developed by Colosi 
and BrowrP^, which is a sum over internal wave modes j 
with horizontal wavenumbers k\. It has the form 



eVi(z,r)= ^2 ^ Vj(z m , h) cos (fai + hr) 



(5) 



i=i 1=1 



where Vj(z;ki) contains the depth dependence and 
weighting of the j th internal wave mode, and there is 
a set of uniformly chosen random numbers {(j>ji) on the 
interval [0; 2tt) that create a random wave field eVi(z, r). 
The upper limits on the summations are set by compu- 
tational, frequency, and environmental considerations^. 
The weighting Vj(z; ki) decays with increasing j and be- 
yond a maximum value makes negligible contributions. 
For the horizontal wavenumber summation over Z, there is 
a minimum and maximum relevant horizontal wavenum- 
ber dependent upon / and thus l max is chosen so that 
there is a sufficiently dense uniformly selected values of 
ki included in the internal wave field construction. See 
Appe ndix [A| for further details, the definition of Vj(z; k{) 
(Eq. ( |A2| )), and the specification of the ki included in the 
summation. 



C. Acoustic modes and unitary propagation 

For a range-independent potential, Vq(z), a general solu- 
tion to the parabolic equation can be written in terms of 
separable solutions 



(*, r; k) = <hn(k)e- ikrE ™rl> m (z-, k) 



(6) 



where t/j rn (z;k) satisfies the Sturm-Liouville eigenvalue 
problem, 



m 1 7,2 



2 dz 2 



(7) 



with ip rn (z;k) and E m the eigenfunctions (modes) and 
eigenvalues, respectively. The boundary condition for 
Eq. ([7| is that the eigenmodes vanish sufficiently rapidly 
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as \z\ oo, the eigenmodes are square integrable, and 
the phase is chosen such that the solutions are real. 
The set of weightings {a m (k)} determine the initial state 
(r = 0). The eigenvalues give each mode a horizontal 
wavenumber kE n for their propagation. Figure [T] illus- 
trates the depth dependence of a few of the first modes. 
Note that slow, adiabatic changes with range can be ac- 
counted for as well in exactly the same manner. There is 
no loss in generality by considering the vertical confine- 
ment to be Vo(z). 

The modes are orthonormal and serve as a complete 
basis for arbitrary wavefields in depth. The scattering 
between modes has been analyzed by looking at both the 
range evolution of the c oupling coefficients 24 and of the 
mode amplitude d 27 * 32 * 3 ^ the former representing scatter- 
ing between mode m and mode n during the propaga- 
tion from internal wave modes. Information about the 
acoustic scattering in each component of the modal basis 
is found by propagating each mode i/j n (z;k) as a con- 
tinuous wave using the full parabolic equation Eq. 0, 
denoted as ^f n (z,r;k), and projecting onto the set of 
ifj m (z;k). Mathematically, the unitary propagation ma- 
trix elements U mn (r; k) are 



^ n (z,r;k)^(z;k)dz . 



(8) 



These elements represent the probability amplitude of 
making a transition from mode n to mode m. Since prop- 
agation with the parabolic equation preserves the norm 
and the orthonormality of two initially orthonormal wave 
fields, this matrix must be a unitary matrix. Therefore 
UU^ = /, where W is the Hermitian conjugate (the com- 
plex conjugate transpose of matrix U) and I is the iden- 
tity matrix. If propagated in the unperturbed system 
with potential Vq(z) (i.e. internal waves are turned off, 
eV\ = 0), each propagated mode \I/ n (z,r; k) accumulates 



only a phase, i.e. ^ n (z,r;k) = e 



-ikrE n 



ip n (z; fc), and the 



unitary matrix is diagonal. If propagated using the inter- 
nal wave perturbation, each propagated mode *ff n (z, r; k) 
mixes into other modes. The amount and nature of the 
mixing is a consequence of the properties of the internal 
wave field. 



E 
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FIG. 1. Eigenmodes, ip n (z; fc), of the waveguide potential, Vb, 
are plotted versus ocean depth, z for n — 0, 1, 10 at their 
eigenvalues E n . The potential Vb is plotted as a dashed line. 



at each source wavenumber^, 



$(*,r;t) 



(z,r;fc)e 



ikco (^t 



'dfc, 
(9) 

where ko = 27r/o/co is the central source wavenumber 
and at is the standard deviation. A source frequency of 
fo = 75 Hz with a standard deviation 18.75 Hz is used 
to model a source with a 3 dB bandwidth of 37.5 Hz^R 

The range-dependent wavefield can be written in 
terms of the propagated modes as \I/(z,r;fc) = 
a n (k)ty n (z, r; k). Thus, the timefront in terms of 
modes is computed as 



$(z,r',i) 



==Y^ / dfc a n (k)U mn (r;k) 
Kcrfr ZZ J 



^m(z;k)e 



m,n 
ikco ( t 



o-fcor 



(10) 



D. Acoustic timefront 

Source signals that are localized in time, as opposed 
to being continuous waves, generate signals at the down- 
range receivers as a series of pulses with time. Consid- 
ering these signals for all depths forms a structure called 
the acoustic timefront. The construction of the time- 
fronts requires a weighted Fourier transfrom of the solu- 
tions at fixed-fc. A Gaussian-shaped, localized pulse in 
time is used here, which results from a Gaussian weight 
in the Fourier transform. More specifically, the complex 
acoustic timefront $ at a range r is constructed from a 
superposition of the continuous wave arrivals \I/(z,r;fc) 



This is the experimentally measurable quantity from 
which to deduce information about the internal waves 
and sound speed profile. Figure [2] shows timefronts 
propagated with and without internal wave scattering. 
Without scattering, it has well defined, stable branches 
throughout its entire duration. With scattering, there 
are branches with fluctuating intensities and small shifts 
in locations in the early arriving part. Most of the power 
is contained in the latest part of the arrival where the 
intensities still fluctuate, but the branches overlap and 
disappear into a complicated finale. The complicated 
structure is deterministic since for a fixed perturbation, 
the resulting propagation is unique, but over the ensem- 
ble of internal wave realizations, it can be analyzed sta- 
tistically. 
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FIG. 2. Color density plots of timefronts calculated using 
the parabolic equation and weighted Fourier transform, which 
illustrate their basic features for propagation of 1000 km. The 
upper figure is without including internal wave scattering, the 
lower figure with internal waves. 



E. Random matrix theory 

Random matrix theory gives a means to understand 
the statistical properties of wave chaotic systems. As 
such, it can be thought of as a branch of statistical me- 
chanics; see Ref. [14] and references therein for a review of 
applications to wave guide problems and quantum trans- 
port in mesoscopic conductors. A complete theory begins 
with a construction of a random matrix ensemble that 
models the chaotic system of interest sufficiently well so 
as to reproduce all of its statistical properties. Ideally, 
that would be followed by analytic calculations as expec- 
tation values over the ensemble which would allow one 
to predict and describe the system's statistical properties 
without even actually having to construct the ensemble 
numerically. Often such calculations are quite sophisti- 
cated and one does not know whether they are tractable 
in advance. Even in their absence, the ensemble may be 
quite useful for modeling purposes and for understanding 
what information exists in the system's statistical proper- 
ties. The o riginal ensembles, introduced by Wigner and 
DysorP^^, are featureless and contain no information 
about the system other than what fundamental symme- 
tries are preserved, i.e. time-reversal symmetry, angular 
momentum!^. For long-range ocean acoustics, some in- 



formation must be built into the ensembles, although one 
does not know, a priori, which information and how to 
incorporate it. That process begins in the next section. 

A critical concept is that random matrix ensembles 
possess a generalized ergodic property. This crudely 
means that in the limit of increasing matrix dimensional- 
ity, each individual system's statistical properties match 
those of the ensemble itself. For the standard ensem- 
bles of Wigner and Dyson, this was proven by Pande}El 
This is a stronger feature than is needed for the ocean 
acoustics system as repeated measurements can be taken 
hours apart by which time the precise state of the ocean 
changes significantly enough so as to be able to gener- 
ate roughly "independent measurements." Nevertheless, 
even a single acoustic field propagation is statistically 
likely to possess the structure and fluctuation properties 
(sampled properly within the timefront) inherent in an 
ensemble of such independent measurements, after care- 
ful accounting for diurnal and seasonal secular behaviors. 

In order to apply random matrix theory fully, it is nec- 
essary to be precise about the statistical measures to be 
calculated. For a timefront, many quantities have been 
discussed such as the wander in the branches from one re- 
alization to the next, the statistical fluctuations in ampli- 
tudes in the finale or somewhere on a particular branch, 
the spreading within a branch, and any bias in the time- 
front to be advanced or retarded relative to the pertur- 
bation free case 2 . It is of interest to understand the fre- 
quency, range, and perturbation strength dependencies of 
the statistical measures. As we are not attempting here 
analytical ensemble calculations, we do not give exact 
statistical measure definitions for this study, but rather 
just recognize that it would be necessary to begin such 
calculations. 

Finally, information is another critical concept. In the 
standard ensembles, only information about fundamen- 
tal symmetries is included. As an example, if the system 
is time reversal invariant, the random Hermit ian matrix 
representing the system Hamiltonian has real Gaussian 
random entries for the matrix elements. For systems in 
which time-reversal is strongly broken, the matrix ele- 
ments are complex Gaussian random entries. This is 
insufficient for ocean acoustics as at least some limited 
information about the nature of the perturbation has to 
be taken into account. A method for constructing ap- 
propriate random matrix ensembles is given in the next 
section. 



III. RANDOM MATRIX ENSEMBLE CONSTRUCTION 

The unitary propagation matrix of Eq. ([8| is the most 
suitable object for which to construct a random matrix 
ensemble. Our construction^, here given in greater de- 
tail, begins with the property of multiplicative unitary 
propagation in which the product of two unitary propa- 
gators gives the propagator for the sum of their ranges. 
The idea is to determine a range of propagation, r^, over 
which the matrix f/j(r^; k) initially gains certain statisti- 
cal properties, specified ahead. The Uj^^k) becomes a 
building block propagatoi^ in which ideally each is sta- 
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tistically independent of its predecessors and successors. 
The block is quite reminiscent of the slices used in the 
Chernov method^. The subscript j indicates then that 
this unitary propagator is one particular realization of the 
building block random ensemble. A single member of the 
full random matrix ensemble is given by the product of 
N statistically, independently-drawn, building blocks as 



N 



U rmt (Nr b ',k) = Y[Uj(r h ;k) , 



(11) 



where the full range of propagation is given by Nr^. It 
remains to be determined what an optimal range for 
is and what are the statistical laws governing the ma- 
trix elements of Uj^; k). The subscript rmt is added to 
distinguish this ensemble member from the unitary prop- 
agator that results from solving the parabolic equation 
for a specific potential realization. In the case that the 
total range of propagation is just r^, then the full random 
matrix ensemble reduces to the individual building block 
random matrix ensemble. 

Ideally, the propagation range is as short as possible 
subject to the constraint that the physical system gener- 
ates {Uj(r^ k)} with the statistical properties: i) correla- 
tions between Ujir^ k) and /7j+ m (r5; k) are small for any 
nonzero integer m, ii) the phases of each matrix element 
[Uj(rb]k)] mn are largely randomized, and hi) dynamical 
correlations between neighboring matrix elements within 
a given Ujir^k) are minimal. All the matrix elements 
of Uj(rb m , k) are correlated due to the constraints of uni- 
tarity. However, the dynamical correlations just men- 
tioned are in addition to those. If is too short in range, 
there is little scattering, ^(r^; k) is nearly diagonal, and 
the matrix elements are highly correlated with seemingly 
little randomness. In a quantum/wave chaotic system, 
as propagation range increases, these statistical proper- 
ties do appear. On the other hand, if is too long, 
one may miss universal properties that emerge from the 
product structure of the propagation. Just as products 
of independent random variables generate lognormally 
distributed quantities due to a central limit theorem (un- 
der certain conditions), products of independently chosen 
structured random matrices may also lead to universal 
laws. In addition, for too long a range, there would be 
no way to attempt using range-dependent perturbation 
theory to guide the ensemble construction. 

For very short propagation ranges, perturbation theory 
gives an approximation for the propagator as 

U(r;k)^A(r;k)^I-iekJ^ dr / V/(r / \k)j , (12) 

where the elements of the diagonal matrix A are given by 



A mn (r;/c) = S mn ex-p(-ikE m r) . 



(13) 



and Vi(r', k) is the operator corresponding to V\{z,r) in 
the interaction picture and whose modal matrix elements 
are given by 



This is similar to the coupling matrix of Ref. 24. A(r; k) is 
the unitary propagator if the perturbation is absent and 
only the wave guide potential taken into account. To 
include adiabatic range dependence would require that 
the argument of the exponential in A mn (r; k) be an inte- 
gral over range as E m would vary slowly and to incorpo- 
rate the r'-dependence of [i£ m , E n , ipm(z; k), ip n (z; k)] in 
I', (/•':/.•!. 

First order perturbation theory does not maintain uni- 
tarity. However, there are multiple ways to create unitary 
propagators, which can be connected to the perturbation 
theory result. One convenient approach is to apply the 
Cayley transformation such that as e — >• 0, the diagonal 
result is recovered. For any Hermit ian operator A and 
identity matrix /, 



U = A{I + ieA)- 1 (I-ieA) , 



(15) 



U is guaranteed to be unitary. The matrix elements of A 
are given in Eq. (13). The key is to Taylor series expand 



the inverse operator and connect the statistical properties 
of the A matrix elements with the statistical properties 
of the matrix elements of Vj by using Eq. ( 12 ). This gives 



A= 2j dr W;*) 



(16) 



if the expansions of U and U (r&; k) are equated. Actually, 
for the construction of the random matrix ensemble, one 



jk(E m -E n y 



'J &zip* m {z\k)ip n {z;k)Vi{z,r') 
(14) 



does not equate A to the right hand side of Eq. (16), one 
equates the probability density of each A matrix element 
in the random matrix ensemble to the probability density 
of a matrix element of the right hand side evaluated over 
the ensemble of internal wave fields. As long as e is small, 
this construction restores unitarity while altering very 
little the precise fluctuation properties of the individual 
matrix elements predicted by the integral over V/(r; fc). 

For the ocean's internal waves, the perturbation is 
much stronger near the surface than deep in the ocean 
to such an extent that an approximation is sometimes 
invoked in which there is an abrupt perturbation near 
the surface and everywhere else it vanishes. In that ap- 
proximation, each cycle of the wave has a single dom- 
inant independent perturbation. This approximation is 
not needed, but it leads to the idea of equating roughly 
with the cycle length of about 50 km. First, perturbation 
theory should handle the perturbation through a cycle 
as that implies a single main interaction of the acous- 
tic and internal waves. As the perturbation has a very 
complicated, seemingly random, appearance, it may gen- 
erate more or less random matrix element phases as de- 
sired in (ii) above. Second, the next cycle should lead to 
a roughly independent interaction as the internal waves 
separated by 50 km are for the most part independent. 
That would satisfy (i) above. Finally, all parts of the 
wave evolution have interacted with the internal waves 
near the surface. This would seem to be essential (neces- 
sary, though not guaranteed to be sufficient) for reducing 
dynamical correlations [(hi) above]. In fact, as shown 
in Refs. [19j |20j this is the least well satisfied criterion. 
'We shall neglect residual dynamical correlations to begin 
with, but realize that if this approximation is not good, 
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we will find statistical deviations ahead. If however, none 
show up, then that indicates that the information con- 
tained in these correlations is effectively being lost in the 
long range propagation. In addition to that knowledge, 
it is not necessary to return to the ensemble construction 
to attempt to incorporate them. 

By propagating the parabolic equation or by evalu- 
ating Eq. ( [16] ), it is found that the off-diagonal matrix 
elements A mn (m ^ n) do, in fact, behave as though 
their phases are random by about 50 km of propagation. 
They also behave like complex Gaussian random num- 
bers. Therefore, only their variance needs to be specified 



phase ex.-p(—ikE m r) that arises from the eigenmodes of 
the unperturbed system. In principle, it has both a small 
drift (from higher order terms in e) and a fluctuating 
part. Here, we incorporate only a fluctuating part in the 
random matrix ensemble. 



IV. ACOUSTIC TIMEFRONTS 



In order to carry out the ensemble construction of time- 



(the mean of Eq. (16) is zero) 



^-ran(^) — &rnn{k)Zrnn{k) 



(17) 



where z mn (k) is a complex Gaussian random variable of 
zero mean and unit variance (with (z mn (k)z rn ' n ' (fc)) = 
$mm'$nri)i an d the /c-dependence has been emphasized 
because it must be accounted for in the construction of 
a timefront. For m = n, the Gaussian random variable, 
z, is real with zero mean and unit variance, but other- 
wise Eq. ( fl7| ) applies. The variance can be calculated 
analytically using Eqs. (14), (16) and the expressions in 
Appendix [Aj which give 



2 2 
mn 



k 2 r 2 ^ 



sine 



sine 2 ' UJm ' nTb 



£K(*)f 



(18) 



where 



otl = k(E m - E n ) ± h 



(19) 



and sinc(x) = sinx/x. The sine- functions can be thought 
of as smeared ^-functions that pick out the most im- 
portant horizontal wave vector components for that par- 
ticular matrix element, which is crudely a function of 

\m — n\. The J2j |^m'n(^)| determines the weighting and 
is calculated by numerically solving for the modal eigen- 
functions. For some potentials, such as the harmonic 
oscillator, the elements could be calculated analytically. 
This expression for the variance generates building-block, 
unitary propagation matrices that have matrix element 
statistics which agree very well with those from solutions 
of the parabolic equation^ 



The integrand of Eq. (14) has three oscillatory compo- 



nents. In order for the integral not to mostly self-cancel, 
there must be an oscillatory component of r) that 
matches to some extent differences between the oscilla- 
tions of ip m (z]k) and ip n (z;k). The higher frequency 
components of Vi(z,r) fall off in intensity, and there- 
fore one anticipates that as the difference \m — n\ grows, 
the matrix elements must decay in strength. This was 
found in Ref. 19 to result in a mostly power-law decay as 
a function of \m — n\ and to be roughly independent of 
the mean index \m + n\/2. 

Perturbatively, the diagonal element A mrn (k) is mostly 
responsible for introducing a dephasing away from the 



fronts according to Eqs. (10), (11), there is one final cor- 
relation to consider. The value of A mn (k) deduced from 
a dynamical system, say using Eq. ( [l6| , at two different 
values of the wave vector k could be correlated. In fact, 
at a minimum there must be very strong correlations for 
two values of k which closely approach each other. This 
alone implies that the values of z mn (ki) and z rnn (k2) for 
a single member of the ensemble (corresponding vaguely 
to a single realization of the internal wave field) should 
be taken as being nearly identical for very small \k\ — &2I 
with some fall off as the difference grows (not account- 
ing for further types of correlations). A priori, the fall 
off could also depend on \m — n\ as well. Numerically, 
there is considerable decay of correlation over the range 
of wave vectors represented in a timefront and there did 
not appear to be a simple rule to describe its behavior. 
Nevertheless, it turns out that the gross approximation 
of taking z mn (k) identical for all values of the wave vec- 
tor is a far better starting point than the opposite "white 
noise" limit in which each z mn {k) is taken as completely 
statistically independent. The results shown below use 
the former, perfectly correlated choice. Again, if statis- 
tical deviations are found that can be tracked back to 
this crude approximation, then that information can be 
deduced from the propagation. On the other hand, if one 
does not see statistical deviations, then the information 
about the correlation decay as a function of \k\ — A^l is 
being lost by the propagation. 



A. Sample timefronts 

We compare the appearance of a single timefront gen- 
erated using the random matrix ensemble technique with 
another generated using the parabolic equation. The 
question is whether the two timefronts would be simi- 
lar in their overall structure, but not in the details of 
their fluctuations. The former (or ensemble calculation) 
proceeds by applying Eqs. ([l0|, ([II]), ([l5]), ( fl7p9l ) follow- 
ing the prescriptions and choices described in the text. 
The latter (or parabolic equation) calculation relies on 
Eqs. 0, ( |3ff5l ), ([9} for a single realization of the internal 
wave field. 

For the propagation distance of a single building block, 
see the comparison given in Fig.[3j Both realizations span 
roughly the same decibel range and possess similar global 
structures. The branches cover the same locations with 
about the correct intensities. The variations of the acous- 
tic power with depth in the finale are quite similar. The 
fine structures that contain the fluctuations are similar in 
character, but differ in the detailed values. The next test 
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FIG. 3. Comparison of sample timefronts from a random 
matrix ensemble construction and parabolic equation propa- 
gation for 50 km. The magnitudes |$| of the timefront are 
plotted in decibels with respect to the largest value of |<3>| as 
a color density plot in depth z and time t. The upper plot is 
a sample |4>| generated from a single member of the ensem- 
ble model. The lower plot is a sample |<3>| from propagation 
through a single internal wave field. 
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FIG. 4. Comparison of sample timefronts from a random 
matrix ensemble construction and parabolic equation propa- 
gation for 250 km. The magnitudes |4>| of the timefront are 
plotted in decibels with respect to the largest value of |4>| as 
a ca color density plot in depth z and time t. The upper plot 
is a sample |$| generated from a single member of the ensem- 
ble model. The lower plot is a sample |<3>| from propagation 
through a single internal wave field is shown. 



of the random matrix ensemble, in particular the build- 
ing block structure, is to compare timefronts for propa- 
gation to longer ranges, which tests Eq. (11); see Fig. [4] 
For this particular case, globally, the structure is again 
similar, however the fluctuations appear to be somewhat 
stronger in the ensemble construction. Whether or not 
this is a part of the statistical variation to be expected 
or an excess in the ensemble construction requires evalu- 
ating ensemble statistical measures. 



B. Average intensity timefronts 

Consider the average intensity 38 (7(z,r;£)) where the 
brackets (...) indicate the ensemble average for the ran- 
dom matrix ensemble and the average over the internal 
wave realizations for the parabolic equation propagation, 



(/(W)) = ^f>*(^;*)l 5 

1=1 



(20) 



Here TV is the number of internal wave realizations or 
random matrix ensemble members and is chosen large 
enough to have small statistical error (N between 1000 
and 2500 in our calculations); the subscript I indicates 
which realization. The average intensity timefronts from 
the random matrix ensemble and parabolic equation 
propagations for 50 km are compared in Fig. [5] In these 
figures, one sees the two previously mentioned global 
structures, the branches and the complicated finale. The 
branches are almost perfectly similar to their appearance 
in a timefront unperturbed by internal waves, other than 
perhaps a small drift in time. The finale is a region 
of overlapping branches, so much so that there is little 
trace of their origins. The average intensities of the ran- 
dom matrix ensemble and parabolic equation are quite 
similar, although there appears to be a few small differ- 
ences in the finale, such as extra mixing might create, i.e. 
more acoustic power at lower depths. Otherwise, there 
is enough similarity between the two constructions that 
the random matrix ensemble could be used to replace 
the parabolic equation as far as the average intensity is 
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FIG. 5. The average intensity timefront (I) = (|$| 2 ) is shown 
as a color density plot with depth z and time t. (Upper) 
The average is taken over timefronts resulting from ensemble 
calculations utilizing 1000 independent random matrices for 
50 km. (Lower) The average is taken over timefronts from the 
propagation through 2500 independent sound speed models to 
50 km. 



concerned. 

A finer look at the comparison is shown in Fig.[6]by fix- 
ing the depth and range, and plotting the resultant time 
traces as a function of time (i.e. section through the time- 
front). The ensemble captures the locations and magni- 
tudes of the branches exceedingly well. In the gaps be- 
tween the branches, there is a background determined by 
the internal wave scattering, i.e. mode-mixing. It peaks 
near (z,t) = ( — l,r/co) and decays in both depth and 
time from this peak. The rate of decay in both depth and 
time of this structure gives a signature for the multiple 
scattering of the internal waves. It appears to be weakly 
enhanced in the random matrix ensemble results rela- 
tive to the parabolic equation propagation results. The 
enhancement appears to increase slowly with increasing 
depth of the trace and possibly with decreasing time. As 
a result, the earlier branches are not quite as well re- 
solved. The origin of this weak distinction between the 
random matrix timefront construction and the parabolic 
equation one has not yet been identified. One possible 
explanation of this weak discrepancy is the approxima- 
tion of perfect coherence as the value of the wave vector 
k changes (z mn (ki) — z rnn (k2) for all fci, A^). This would 
exaggerate the effect of the perturbation. Less likely, but 
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FIG. 6. Traces of the average intensity are shown for several 
final depths z (0.5, 1.0, 1.5, 2.0, 2.5, 3.0 km, respectively) 
as a function of time t. The average of the ensemble model 
timefronts (solid lines) is taken from 1000 independent ran- 
dom matrix ensemble members for 50 km. The average of the 
propagated timefronts (dotted lines) is taken from the prop- 
agation through 2500 independent sound speed models to 50 
km. 



also a possible explanation, is ignoring the weak dynam- 
ical correlations that must be present in any continuous 
dynamical system. Identifying the precise cause of the 
small differences and improving the ensemble timefront 
construction is beyond the scope of this paper and is 
left for future work. Regardless of the origin, it may be 
possible to make a small rescaling effectively reducing 
the {cr^ n (&)} from their perturbation theory values and 
maintain just as simple an ensemble as the current one. 

To test how the multiplicative structure of the ensem- 
ble fares for the average intensity, one must again look at 
longer ranges. In Fig. [7] an average timefront comparison 
is shown for propagation to 1000 km. The average in- 
tensity timefronts from the ensemble model capture the 
finale region very well both in decay with time and with 
depth, and the positions of the branches are excellent. 
However, the excess scattering results in the branches 
being less visible in the color density plot. Looking at 
the time traces in Fig. [8] gives a sharper view just as it 
did for Fig. [6j In the upper plot, the internal waves are 
turned off to illustrate that in between the branches, the 
signal is many orders of magnitude decreased and one ef- 
fect of internal waves is to infill a structure between the 
branches. In the lower plot, the random matrix ensemble 
captures the fill-in very well, but shows some excess in- 
tensity in the gaps relative to the parabolic equation and 
the magnitudes of the branches are reduced somewhat 
though still in the correct positions. We suspect that im- 
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FIG. 7. The average intensity timefront (I) = (|<£| 2 ) is shown 
as a color density plot with depth z and time t. (Upper) 
The average intensity timefront results from 25 independent 
random matrix ensemble members propagated for 1000 km. 
(Lower) The average intensity timefront is from the propaga- 
tion through 28 independent sound speed realizations to 1000 
km. 



proving the building block ensemble would alleviate this 
slight discrepancy at longer ranges and it is not due to 
the the complete independence assumption of adjacent 
building blocks. 



C. The mixing front 

In the ensemble averaged intensity, (I(z,r;t)}, the 
branches are largely unchanged from the unperturbed 
system. That suggests looking at the ensemble averaged 
intensity change 



5I e (z, r; t) = (I e (z, r; t)) - I e =o(z, r; t) 



(21) 



whose functional dependence reflects the scattering due 
to the internal waves. This is related to the mode-mixing, 
albeit over a range of wave vectors, because the early ar- 
rivals represent wave energy in the higher modes (steeper 
propagation - reverse dispersion case). If there were com- 
plete mixing the branches would disappear, and the early 
and late parts of the timefronts would be statistically and 
structurally indistinguishable. As this is not the case, it 
is possible to deduce information about the mode-mixing 
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FIG. 8. Traces of the average intensity (I) = (|0| 2 ) are shown 
for a single final depth z — 1.5 km as a function of time t. 
(Upper) A comparison of the average propagated timefront 
traces with and without internal waves (solid line without 
internal waves). (Lower) The average of the ensemble time- 
fronts (solid line) is from 1000 independent members for 1000 
km. The average of the propagated timefronts (dotted line) is 
from 2500 independent internal wave realizations propagated 
to 1000 km. 



from 5I e (z, r;t). We call this the mixing front, whose de- 
cay toward the early part of the arrival (as a function 
of r/co — t) reflects the modal decomposition of the ini- 
tial signal convolved with the mixing from the scatter- 
ing. For short enough propagation ranges, this should 
be dominated by the modal decomposition of the initial 
signal, but as the range increases the power-law decay of 
the unitary matrix elements as a function of \m — n\ must 
eventually dominate. More generally, the average inten- 
sity of a large number of well-separated-in-time, long- 
range, ocean acoustic measurements would allow one to 
deduce information about the sum of all the scattering 
processes relevant to those measurements, i.e. that is if 
the experimental noise levels are sufficiently controlled 
that the mixing front can be made to emerge from the 
average. 

The random matrix ensemble for the acoustic propa- 
gation can be used to derive an expression for 5I e (z, r; t) 
without having to perform full propagations over a large 
ensemble of internal wave fields. Starting by approximat- 



Hegewisch et al.: Constructing acoustic timefronts... 9 



32.8 



33.3 
t (sec) 



33.8 



-30 



-20 -10 
ln(<l> ) 







FIG. 9. The solid line depicts the mixing front 8I £ (z,r;t). 
The peaks have disappeared. The average intensity without 
internal waves shows where the peaks of the full timefront 
would be. 



ing the unitary propagation matrix elements as 



ikr b Err 



(Sr, 



(22) 



and performing the ensemble averages using the statis- 
tical properties of z min , the ensemble averaged intensity 
change, i.e. the 'mixing front', takes the form 



SL 



+oo 



dk e-M*-^) 



exp 



~(k - kg) 2 

2*2 



FIG. 10. The average intensity (/) of the acoustic timefront 
is plotted with depth z at the time t=r/co for the ensem- 
ble model prediction (I e ) from Eq. pi] ) (Solid line), for a 
statistical average of intensities from several perturbed prop- 
agations (Dotted line) and for the unperturbed propagation 
ie=o (Dashed line). 



taken into account. This prediction closely approximates 
the average intensity computed from perturbed propaga- 
tions. Thus the random matrix ensemble also gives an 
explanation of the increase in vertical extent of the ar- 
rivals of finales of the perturbed timefronts over those of 
the unperturbed timefronts. 



2a Am Jk)e-* krbE ™a n (k)ip m (z;k)f 



(23) 



The mixing front for a single depth is illustrated in Fig. [9] 
The effect of the mixing front on the timefront is to hide 
most of the structure of the unperturbed arrivals in the 
finale region, with only a slight broadening of the unper- 
turbed early arrivals. This creates the appearance of a 
smearing of the arrivals in the finale region. Thus, the 
random matrix ensemble model gives an explanation of 
the appearance of the arrivals smearing in the finale re- 
gion and why the unperturbed arrivals are still visible in 
the early arrivals of the timefront. As a practical consid- 
eration, experimental determination of the mixing front 
would be enhanced by a better signal-to-noise for shorter 
range propagations. However, longer range propagation 
would be more dominated by the power-law behaviors 
and the mixing front would emerge more clearly. A bal- 
ance of these competing factors would determine the best 
ranges for testing mixing front predictions with measure- 
ments. 



Figure [10| illustrates the depth dependence of the mix- 
ing front of Eq. (21). As seen in the figure, SI e results 



in a greater intensity for the arrivals in the finale region 
(compared to the unperturbed propagation) and a slower 
decay in depth. The decay in depth may transition from 
an exponential to a power-law decay as internal waves are 



V. CONCLUSIONS 

The random matrix ensemble for the unitary propa- 
gation matrices capture the statistical properties of the 
propagation^ quite well without fitting parameters. The 
ensemble is derived using first order wave perturbation 
theory. The ensemble predicts the presence of a mixing 
front in the acoustic timefronts. This front is consistent 
with the long-known presence of a smeared region of ar- 
rivals in the finale region of the timefronts and the persis- 
tence of branches in the early timefront. Its characteristic 
features are the decay with depth and with time from a 
peak near the horizontal arrival. These features seen in 
the ensemble agree with parabolic equation propagation. 

Though perturbation theory alone would only be able 
to make predictions about the unitary propagation ma- 
trices to very short ranges, the statistical properties pre- 
dicted using perturbation theory seem to be more robust 
over the course of the propagation to much longer ranges 
(i.e. 3000 km). The successful use of first order pertur- 
bation theory in capturing useful statistics for the long 
range propagation suggest that perturbation theory may 
have statistical use for other fields as well. 

The realization that internal wave effects cause a mix- 
ing front in the acoustic timefront of arrivals may lead 
to a new way of deducing information about the internal 
state of the ocean. This front has characteristic features 
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in depth and time that depend on the properties of the 
mode mixing. Their parameters are controlled by the 
strength and structure of the internal wave sound speed 
perturbations. For example, measurements of the ver- 
tical extent of arrivals in the finale region have already 
been analyzecP^Hl. Comparing to the predicted proper- 
ties of a mixing front may give greater understanding of 
the ocean environment. 

There are many possible directions for future research. 
More careful identification of information lost and re- 
tained and what effects remain would be very interesting. 
There is also the possibility of extending the construction 
of random matrix ensembles to broader classes of ocean 
acoustics systems where absorption plays a role. More 
work should also be done with the fluctuation properties 
of the timefronts and attempts made to have analytic 
random matrix predictions for their behaviors. 
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APPENDIX A: IMPLEMENTATION OF INTERNAL 
WAVE MODEL 

The efficient numerical scheme devised by Colosi and 
Brown® generates a statistical ensemble of internal wave 
effects for the sound speed model through the equation: 



eVi(*,r) = 2 ^ 2 ^7V 2 A /^exp(-3,/25) (Al) 



g 7T 



M 



jmax In 



b max / j- 

X J2 Yl sin (M(*0) \ ■ ;2+ fc '-2 C0S + fc ' r ) ' 
J I J* 



3 = 1 1 = 1 



where the depth dependence and weighting of the j th 
internal wave mode is 



9 7T 
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1 , 1 1 + i + l 



/? 2 + i 2 + i)§ y^j^i-i) ' 

(A2) 

and £(z) = e~ z l B — e~ H l B . A single random seed gen- 
erates a sequence for the random phases, <j)ji G [0,27r), 
of each internal wave with vertical mode, j, and hori- 
zontal wavenumber, k\. The ensemble has the statisti- 
cal properties of the Garrett-Munk spectrum^. The full 
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